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Abstract. We consider pairing in a two-component atomic Fermi gas, in a three- 
dimensional optical lattice, when the components have unequal densities, i.e. the gas 
is polarized. We show that a superfluid where the translational symmetry is broken by a 
finite Cooper pair momentum, namely an FFLO-type state, minimizes the Helmholtz 
free energy of the system. We demonstrate that such a state is clearly visible in 
the observable momentum distribution of the atoms, and analyze the dependence of 
the order parameter and the momentum distribution on the filling fraction and the 
interaction strength. 
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1. Introduction 

Newly realized strongly interacting Fermi gases [HEllSllllEllSllZllHlini have enabled 
an intriguing line of research into the nature of fermion superfluidity. In ultra-cold 
degenerate gases the system is composed of atoms in different internal states. The 
number of atoms in different states can be experimentally controlled and this has enabled 
first experimental studies of polarized trapped Fermi gases |TDJ |TTJ EE2] Imbalanced 
atomic gases have recently inspired a number of theoretical works [T3| ITU ITK| ITo^ IT7| 

immi23injl22l23ElE21l2n],in addition to the wide literature on fermion pairing 
with unequal chemical potentials or number densities in the fields of condensed matter, 
nuclear and high-energy physics (for a review of some of these studies, see |27|). 

While several parameters of a harmonically trapped Fermi gas can be controlled 
experimentally, the inhomogeneity is often crucial and some parameters, such as the 
atomic mass, cannot be changed. The use of optical lattices provides unprecedented 
tunability and, for example, the effective masses of atoms can be changed by simply 
changing the laser intensities. Furthermore, in the center of the optical lattice the 
effects due to the harmonic trapping can be weak. This versatile tool has been recently 
widely used and has enabled, among other things, the observation of the superfluid-Mott 
insulator quantum phase transition [2E| i n a cloud a bosons. Use of optical lattices is not 
restricted by quantum statistics and indeed fermions have been studied experimentally 
in one-dimensional |2TH l3~U] as well as in three-dimensional |3~Tf I32[ 13*3*] optical lattices. In 
addition, the first experimental results of Bose-Fermi mixtures in optical lattices were 
recently reported |34j . 

In this paper we investigate the properties of polarized Fermi gases in a three 
dimensional optical lattice. In particular, we focus on the states of constant density and 
compare in detail the energetics of the breached pair (BP) states (33 EH1 E3 EH] and 
the simplest variant of the FFLO states p""")} "4T"] in a three-dimensional optical lattice. 
We find that, while the non-zero gap BP state can be a minimum of the Helmholtz free 
energy it can lower its energy by forming pairs with non-zero momentum. 

We show that the lattice dispersion causes a different dependence of the superfluid 
gap on the polarization compared to homogeneous systems. We find that, while the 
non-zero gap BP state can be a minimum of the Helmholtz free energy it can lower 
its energy by breaking the translational symmetry and forming pairs with non-zero 
momentum as in FFLO states. Section 2. reviews the theory used in the numerical 
calculations. In section 3. we discuss the free energy analysis and explain, illustrated 
by numerical examples, why it is necessary to use the Helmholtz free energy in case 
of a system isolated from the environment with respect to particle exchange, as is the 
case with trapped atoms. In section 4. we present the results on the FFLO states: 
we analyze the system behaviour as a function of the lattice filling fraction and the 
interaction strength between the two pairing components. Furthermore, we present 
the momentum distribution of the atoms, showing that this easily observable quantity 
carries a clear signature of the FFLO state. We conclude by a discussion in section 5. 
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2. Theory 

The microscopic Hamiltonian for fermions of two (pseudo) spins f and I, in an external 
potential V(x), is 

H = £ J ^(x) (-|^V 2 + V(x)) M*) d3 *+^ 2 ^ I V? WVl(x)Vi(x)V^(x) rf 3 x, (1) 

where a is the s-wave scattering length. With a sufficiently deep periodic potential, the 
atoms are localized in the minima of the potential, and the system can be described by 
the Fermi-Hubbard Hamiltonian H"T1. 



H - /i T iV T - fl l N l = - (W{nCtn + mdl^n) + £/ E c\ n c\ n C lrl C U 




E +Jm E (4mC T n + c\ m C ln ) . (2) 
( n > m >2/ <n,m) z / 

Here (n, m) x means a nearest neighbour pair in the x-direction. In mean-field theory 
the interaction term in the Hamiltonian is approximated with 

U E a Tn4n a inC Tn = U E ((^In) d ln^n + cln^n ( d ln^n) ~ (c\ n c\ n ) (q n C Tn )) , (3) 
n n 

where the Hartree and Fock terms have been dropped since the former are effectively 
included in the chemical potentials and the latter do not contribute. A general ansatz 
U (c| n C| n ) = Ae 2iq ' n describes several apparently different states. For equal (pseudo)spin 
number densities, q = gives the standard Bardeen-Cooper-Schrieffer (BCS) solution. 
For unequal densities q = gives the Breached Pairing (BP) solution and non-zero 
values for q describe a Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state where Cooper 
pairs have a finite momentum, q. Finally, A = naturally describes a state with no 
superfluid, i.e. the normal state. 

Including the chemical potentials in the Hamiltonian it becomes, in momentum 
space, 



IA 



-) , (4) 



H = E (( e k - 4k a Tk + (£k - fii) 4kC|k + Ac{ q+k c{ q _ k + h.c 

k \ / 

where M is the number of lattice sites, /if and /i| are the chemical potentials of the 
different spin species and the lattice dispersion is given by 

e k = 2,4(1 - cos(M)) + 2J y (l - cos(k y d)) + 2J Z (1 - cos(k z d)). (5) 

Here d is the lattice parameter, i.e. the distance between two neighbouring lattice points. 
The parameter U describes the energy associated with the interaction of the particles, 
with negative values corresponding to an attractive interaction. The hopping parameter, 
J is the energy gain corresponding to tunneling between nearest neighbour sites. In a 
lattice, J is essentially the band width. For a more detailed discussion on the parameters 
J and U, see e.g. jl2j. Using the fermionic anticommutators, the Hamiltonian can be 
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H =jr I Y,\ ( £ k+q - ft) 4k+q £ Tk+q ~ (f-k+q ~ Hi) (c^fc+qcj-k+q ~ l ) + 
k \ 

c c|q-k + A*C iq _ k C Tq+k - 



Acjq+k' |q- 



This can be written in matrix form as 



"77 



(6) 



k 



^Tk+q ^l-k+q 



Ck+q - 

A* 



A 

(e-k+q - Hi) 



C Tq+k 
C |q-k 



(7) 



Because A was chosen as the amplitude of the order parameter, it is a real number, 
which simplifies the expressions. The second sum in (J7J) is just a constant, but it is 
important for the calculation of free energies. The eigenvalues of the matrix part (i.e. 
the quasiparticle energies) are 



E: 



Hi — Hi . e k+q - C-k+q 



Ck+q + e_k+ q Hi + Hi 



M, ± = — j -V ^ j ■ —J + A'.(8) 

With the introduction of the single particle energies = e k — Hi an d £|k = e k — Hit 
this can be written in a simpler form 



^k,q,± = = ± 



\ 



£?k+q + ^i-k+q 



A 2 



(9) 



Unequal densities introduce unequal chemical potentials, which destroys the particle- 
hole symmetry that exists in the standard BCS theory. A suitable Bogoliubov 
transformation, 



/ c Tq+k \ _ I 

{ 4q-k ) ~ [ 

with the coefficients given by 



u 



k,q 



k,q 



ill 



M k ,q ^k,q 
-fk,q Mk,q 



7Tk,q 



(10) 



a+A 2 



A 



Mk,q^k,q 



2,/ei + A 2 



where £a = (£tk+q + £l-k+q)/2, diagonalizes the Hamiltonian to 



H 



M 



E 



k L 



7Tk,q 7i-k,q 



-^k,q,H 






"^k,q - 



7Tk,q 
7i-k,q 



£|-k+q 



A 2 
17 



(11) 



(12) 



Here the 7 CT ,k,q are the quasiparticle operators and they obey the Fermi statistics. 
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2.1. Self- consistent equations for A, and 

The ansatz U (q n cj n ) = Ae 24q n implies A = U ^ k (c|_k+qCtk+q) /M, where M is the 
total number of lattice sites. With the Bogoliubov transformation ([ID}) , this can be 
written as 

A = ^ £ (q-k+qC T k+q) = ^S«k,qVk,q(/(^k,q,-) ~ /(^k,q,+ )) 
k k 

^ ^ /(-^k.q ,-) ~ /(^k,q, + ) / 13 \ 

M k 270T^ 



s(A) = A- V 7 l " k ' q '-/ ~ 7 - A = 0. (16) 



The following equations hold for the particle numbers: 

= E (4 £ Tk) = E <q/(£k,q, + ) + <q/(^k,q,-) (14) 
k k 

Nl = E(4 £ ik) = E<q/(-^, q) -)+< q /(-^k,q J+ ). (15) 
k k 

Equations (|T3|) . (jEJ), and (|TK|) are self- consistent equations in the sense that they are 
always satisfied by the stable macroscopic state. Especially the gap equation, (JT3J), can 
have many solutions, of which the one with the lowest free energy is stable. For example, 
for equal populations of the two components, even below the critical temperature, the 
gap equation always has the trivial A = solution, in addition to the standard BCS 
solution. However, the BCS solution has lower free energy and is therefore stable. 

It is useful to define a gap function g(A) that vanishes at the correct value of A. 
This reduces dTJ) to g(A) = 0. 

U y /(£k,q,-) ~ /(g k,q, + ) 

In optical lattices, the relevant quantity, instead of the total number of particles in 
one component, is the filling fraction of one component. The filling fraction is defined 
as the number of particles divided by the number of lattice sites. We denote filling 
fractions as /f = N^/M and fi = NJM. Because of the Pauli exclusion principle, only 
one fermion of each kind can fit on the lowest Hubbard band at each site, thereby having 
both /f and /j equal to 1 means having a full lattice, which is a band insulator. Putting 
three atoms in the same site or on the same momentum state would mean populating 
higher energy bands, which involves a significant cost in energy. For superfluidity, the 
optimum setting, giving the highest value for A, is /j = 0.5 and f\ = 0.5|13j. For 
population imbalance, it is useful to define the polarization: 

p = f -frr- (it) 

h + h 

2.2. Experimental parameters 

In our numerical calculations we consider 6 Li atoms trapped in a lattice created by 
A = 1030 nm lasers. The wavelength is related to the lattice parameter d by d = A/2. 
The lattice height is 2.5Er, where the recoil energy is h 2 (2n / A) 2 /2m. In the present 
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paper calculations are done in zero temperature. Except where explicitly noted, we use 
— 1000 Bohr radii for the scattering length a. Typical lattice sizes we have used are 
64 x 64 x 64 and 128 x 128 x 128. We calculate the finite lattice sums explicitly, which 
limits all the vectors in our numerical analysis to a set of discrete lattice points. 



3. Free energy analysis 



The relevant free energy of a Fermi gas depends on the physical system in question. 
We start by considering the grand potential Q = — 4 InZc, where Zq = Tre _/3H , is the 
partition function of the grand canonical ensemble. The partition function is 



2 G = Tre^ = n 1 



^k,q, 



1 + e* 



^k.q,- 



-I3C 



where C corresponds to the constant part of the Hamiltonian and j3 = 1/k^T. The 
grand potential is 



n 



1 

M 



E 



k+q 



E: 



k,q, 



A2 

U 



In l + e" 



^k,q, 



+ In 1 + & 



"k.q, 



At low temperatures this is independent of (3 and becomes 



n 



1 

M 



k+q 



E 



k,q, 



U 



+ -Ek,q, + @( — -Ek,q, + ) — -^k,q -©(-^k.q, 



.(19) 



(20) 



There are two different schemes for treating population imbalanced Fermi gases: 
fixing the particle numbers or fixing the chemical potentials; this choice depends on the 
physical system in question. We are interested in the former, but first we will briefly 
discuss the latter, to enlighten the differences between the two cases. 



3.1. Fixed chemical potentials, q = 

With fixed chemical potentials, the relevant thermodynamic free energy is the grand 
potential fi(A, /if, /if). The extrema of Q, dfl/dA = 0, correspond to solutions of the 
gap equation, JTBJ), i.e. g(A) = 0. Figure [T] shows a typical scenario with the chemical 
potentials fixed at /if = 0.394 and /if = 0.304. The figure shows how the extrema of Q 
coincide with the zeros of g(A). In this scheme, the BP state is unstable because it is a 
local maximum of Q(A). However, both the polarization and the total number of atoms 
change as A changes, as can be seen from figure El Therefore the extrema correspond 
to situations with different total numbers of particles. If the numbers of atoms are to 
stay fixed in the system, this comparison is not valid. 



3. 2. Fixed numbers of particles, q = 

When numbers of particles, instead of the chemical potentials, are fixed, the relevant 
thermodynamical quantity is the Helmholtz free energy, 

F(A, N h ATf) = Q + /if-t + /if -J- = n + /if/f + /if/f . (21) 
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Figure 1. The grand potential per lattice site (green), in the units of recoil energy, 
and g(A)/A (blue) in arbitrary units. The maximum of £1 at A = 0.04 corresponds to 
the BP state with polarization 0.12, and the minimum at A = 0.068 corresponds to the 
BCS state with a zero polarization. The chemical potentials are fixed at /ij = 0.394 
and hi = 0.304. 
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Figure 2. The total filling fraction / = /f + /j, as a function of A, with the chemical 
potentials fixed at ^ = 0.394 and hi — 0.304. 



Now the physical solutions are the minima of F(A, iVj, iVj: dF/dA = gives the 
extrema and d 2 F/dA 2 > defines the minima. These solutions again coincide with 
the zeros of g{A). Moreover, the extrema given by F are the same as given by Q, i.e. 
dQ/dA = dF/dA = 0, but determining which of the extrema are minima, and thereby 
the physical solutions, depends on whether F or Q is used. The solutions where q is 
fixed at zero are the uniform solutions. If the densities of the different components are 
the same, the solution is known as the BCS state and if the densities differ, the solution 
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Figure 3. The Helmholtz free energy F(A) per lattice site (green), in the units of 
recoil energy, and g(A)/A (blue) in arbitrary units. The filling fractions are fixed at 
/l = 0.44 and /j_ = 0.36, so the polarization is 0.1. The minimum of F at A ps 0.045, 
coinciding with the zero of <?(A), corresponds to the BP state. Also shown is 17(A) (red) 
with chemical potentials fixed so that the filling fractions have the values mentioned 
above at the point where g(A) = 0. 



is known as the BP state. With fixed densities, both the BCS (JVj = iVj_) and the 
BP (iVf iVj_) state correspond to a minimum of F(A). We consider there the case 
N^Ni, that is, the BP state. 

Figure El shows the Helmholtz free energy and the gap function with fixed numbers 
of particles (with iVj 7^ iVj,) and q fixed to 0. It is clear that F(A) is minimized with 
a finite A, i.e. the BP state is stable in this consideration, however one has to consider 
also the possibility of a non-zero q. 

4. FFLO-states 

Stability analysis is always limited to some set of states. Many states have been studied 
in isotropic systems. [S3 E01 ESI iSl EI A Monte Carlo study 148] suggests the FFLO 
state to be the ground state in a two-dimensional lattice in the weakly interacting regime. 
Here we study BCS, BP, and single mode FFLO states in three-dimensional lattices. 

When the particle numbers are fixed and the momentum of the Cooper pairs, 2q, 
is allowed to get non-zero values, the translational symmetry is broken and the state 
is FFLO-like. The stable state is now found by minimizing F with respect to A and 
q. Figure |U shows the momentum distributions of the different components along the 
k z = plane in BP and FFLO states. Note that there is a vacant region, or breach, in 
the momentum distribution of the minor component in the BP phase. 
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Figure 4. The momentum distributions /j and /i, on the k z = plane, of the BP 
(top row) and FFLO (bottom row) states. Both states demonstrate a depairing region 
around the Fermi surface, but whereas the region is symmetric in the BP state, it is 
asymmetric in FFLO. The background scattering length is —1500 ao, the total filling 
fraction is 0.1 and polarization 0.2. In this situation, the FFLO state is energetically 
favorable. 

Figure El shows a typical free energy landscape for the FFLO state. The filling 
fractions are fixed at /j- = 0.24 and f\ = 0.16. The figure shows that the minimum 
energy is found with a non-zero A and non-zero q. The BP state is a saddle point and 
the BCS solution is absent since it does not support polarization in zero temperature. 

Figures HI - ^2 show the energy gap A and the FFLO momentum q as a function of 
polarization for different interaction strengths. With increasing polarization, the energy 
gap A decreases and the magnitude of q increases. Whether a critical polarization 
P c , where A vanishes, exists, depends on the total density /j + f± and the interaction 
strength between the atoms, characterized by the scattering length a. With a = —1000 
Bohr radii and total filling fractions between 0.4 and 1.0, the critical polarizations are 
around 0.3. Raising the scattering length to a = —1500 a gives a Pc around 0.6 with 
the same densities. With a = —2000 a , P c is more than 0.9, as can be seen in figure 
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Figure 5. The Helmholtz free energy F(A, q) per lattice site. Here q is in the in- 
direction, i.e. q = (q,0, 0). The filling fractions are fixed at /f = 0.24 and f± = 0.16, 
so that P — 0.2. The units of q are selected so that (32, 32, 32) would correspond to 
the corner of the first Brillouin zone, the R point. 




Figure 6. The order parameter A (blue) and the magnitude of q (red) as a function of 
polarization P. Here the total filling fraction is 0.1, the scattering length is — 1000ao, 
and the lattice size is 128 x 128 x 128. The q is given in reciprocal lattice indices, so 
that 64 on the j/-axis would correspond to the edge of the first Brillouin zone. 




EU 

Figureslol-fTUlshow the q = 0, i.e. the BP phase, with small polarizations. However, 
the discrete steps in the values of q are due to the finite size of the lattice and using 
a larger lattice allows for smaller steps with shorter intervals. We expect these steps, 
and the kinks in the energy gap, to vanish in the limit of larger lattices. The figures 
also show that stronger interactions as well as larger densities lead to higher critical 
polarizations, as is expected. 

Because the Cooper pairs each carry momentum 2q, it would seem that the system 
has total momentum. However, it has been shown for homogenous systems that the 
net momentum is canceled by momentum distributions of the individual components 
canceling the effect of the momentum q This is consistent with figures [T21 - IT71 where 
the total momentum distribution (= /f + /j.) is biased to the direction opposite to 
where q is located. We have also numerically checked that the net momentum is zero. 

Figures ^] - El show the momentum distributions in different states, integrated 
over the k z . The background scattering length is —1000 ao in each figure. The figures 
show clearly the effect of filling fraction, with low filling fractions giving a spherical 
Fermi surface, but higher filling fractions showing deformations caused by the lattice. 
The BP state is visible in the difference of the momentum distributions, — njjc, as a 
symmetric depairing region, whereas the FFLO phase has asymmetric depairing. Note 
that the vacant region in the momentum distribution in the BP phase is present, but 
not visible due to the column integration and also due to the small polarization. 

The FFLO state has been suggested to be observable in the correlations in the atom 
shot noise |5()j. however, as figures IT^I - IT71 show, the FFLO state is already reflected in 
the momentum distributions. 
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Figure 8. The energy gap A and the magnitude of q as a function of polarization 
with the total filling fraction of / = 1.0 and scattering length a = — lOOOao. 




Figure 9. The energy gap A and the magnitude of q as a function of polarization 
with the total filling fraction of / = 0.1 and scattering length a = — 1500ao. 



5. Conclusions 



We have shown that, in the case of interacting fermionic atoms trapped in optical 
lattices, a stable FFLO state can be found. The stability of various non-BCS superfluids 
such as BP state/Sarma state has been a intriguing topic of study and it is known that 
the stability depends for instance on whether the particle number in the system is fixed 
or not. Moreover, one has to always consider also the possibility of a state beyond 
the BP/Sarma state, such as FFLO-type states associated with a non-zero Cooper pair 
momentum. FFLO state has been predicted to lower the system free energy compared to 
the BP state in many cases, however, it is often predicted to occur only in some narrow 
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parameter window. I^Uj According to the analyses we present here, in atomic gases 
trapped in optical lattices, the FFLO state minimizes the relevant free energy, and does 
it for a considerable range of parameters. This is influenced by the fact that the particle 
numbers, not the chemical potentials, are fixed in trapped atomic gases. Furthermore, 
considerable critical polarizations can be achieved, which may be partly aided by lattice 
features. For low filling fractions, e.g. / = 0.1, critical polarizations are P c = 0.075, 
P c = 0.5, and P c = 0.9 for the scattering lengths — 1000a , — 1500a and — 2000a , 
respectively, showing strong dependence on the interactions. Also raising the filling 
fraction increases the critical polarization, for instance for — lOOOao scattering length, 
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Figure 12. The momentum distributions /j + /j,, and /j — integrated over the 
z-direction. Here the total filling fraction is / = 0.1 and P = 0.02. Here q = and the 
momentum difference shows the BP depairing symmetrically along the Fermi surface. 




Figure 13. The momentum distributions f-\ + /j., and /j — /i, integrated over the 
z-direction. Here the total filling fraction is / = 0.1 and P — 0.04. The state is of 
the FFLO type with a finite q = 7t/(128g?), which shows clearly in the momentum 
difference. The depairing region is similar to what has been predicted for homogenous 
systems, see 



P c = 0.3 for / = 0.4 (compared to P c = 0.075 for / = 0.1). Interestingly, however, 
further increase of / does not necessarily make P c to grow, for example P c = 0.3 also for 
/ = 1. This dependence of the critical polarization on the filling fraction requires further 
study and may be useful in comparison of theory to experiments. Stronger interactions 
may bring in additional interesting features, but then one has to consider the validity 
of the single band approximation. We have also showed that the FFLO state is clearly 
reflected in the directly observable momentum distribution of the atoms. 
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Figure 14. The momentum distributions /j + /j,, and /f — integrated over the 
z-direction. Here the total filling fraction is / = 0.4 and P = 0.03. Here again q = 
and the momentum difference shows the BP depairing symmetrically along the Fermi 
surface. 




-1 -1/2 q 1/2 1 -1 -1/2 q 1/2 1 



Figure 15. The momentum distributions /j + and /j — integrated over the 
^-direction. Here the total filling fraction is / = 0.4 and P — 0.15. The state is FFLO 
with q = 27r/(64d) in the ^-direction. 

An important topic of further study is to evaluate the effect of the residual 
harmonic trapping potential which is always associated with optical lattices, despite 
their approximate periodic translational invariance. Issues related to potential phase 
separation have to be clarified. Based on previous experiments on bosonic and fermionic 
atoms in optical lattices, it seems possible to limit the effect of the harmonic potential 
in such a way that essential predictions of a homogeneous lattice model can be observed. 
On the other hand, for strong harmonic trapping, interesting combined lattice and phase 
separation effects may be seen. 
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Figure 16. The momentum distributions /j + f±, and /j — integrated over the 
z-direction. Here the total filling fraction is / = 1.0 and P = 0.01. Here the state is 
BP. 




-1 -1/2 q 1/2 1 -1 -1/2 q 1/2 1 



Figure 17. The momentum distributions f-\ + f±, and /f — /i, integrated over the 
^-direction. Here the total filling fraction is / = 1.0 and P — 0.08. Here the state is 
FFLO with q = 27r/(64d) in the x-direction. 
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